rm(list = ls())

# Install required packages if not already installed, then load
required_packages <- c("haven", "raster", "sp")
new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if (length(new_packages)) install.packages(new_packages, repos = "http://cran.us.r-project.org")

library(haven)   # read_stata
library(raster)  # shapefile
library(sp)      # spatial classes

# Set paths relative to the Reproducibility Package root
prgms_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(prgms_dir)
setwd("..")
work_dir  <- getwd()
rawdata_dir   <- file.path(work_dir, "02_RawData")
finaldata_dir <- file.path(work_dir, "03_FinalData")
figures_dir   <- file.path(work_dir, "05_Figures")

bfa.shp.dir <- file.path(rawdata_dir, "bfa_adm_igb_20200323_shp")
eth.shp.dir <- file.path(rawdata_dir, "eth_adm_csa_bofedb_2021_shp")
mwi.shp.dir <- file.path(rawdata_dir, "mwi_adm_nso_hotosm_20230405_shp")
nga.shp.dir <- file.path(rawdata_dir, "nga_adm_osgof_20190417")
country.list = c("bfa", "eth", "mwi", "nga")
countries <- c("Burkina Faso", "Ethiopia", "Malawi", "Nigeria")
path.list = c(paste(bfa.shp.dir, "/bfa_admbnda_adm0_igb_20200323.shp", sep = "")
              , paste(eth.shp.dir, "/eth_admbnda_adm0_csa_bofedb_itos_2021.shp", sep = "")
              , paste(mwi.shp.dir, "/mwi_admbnda_adm0_nso_hotosm_20230405.shp", sep = "")
              , paste(nga.shp.dir, "/nga_admbnda_adm0_osgof_20190417.shp", sep = "")
)
gps_coord = c()
country_admin0_shapes  = list()
for(i in 1:4){
  country_admin0_shapes[[i]] <- shapefile(path.list[i])
  # gps_<country>.dta files are built by 01_Prgms/gps_coord.do and saved to 03_FinalData/
  gps_coord = rbind(gps_coord, read_stata(file.path(finaldata_dir, paste0("gps_", country.list[i], ".dta"))))
  
  
}
# Set up the plot window
cex = 1.25

png(file.path(figures_dir, "FigureA1.png"), width = 4000, height = 3000, res = 300)
par(mfrow = c(1, 4), mar = c(1, 1, 1, 1))
# Plot each country
for (i in 1:length(country_admin0_shapes)) {
  plot(country_admin0_shapes[[i]], main = countries[i], cex.main = 2)
  points(gps_coord$lon[gps_coord$cc ==i & gps_coord$pnl_urban == 1], 
         gps_coord$lat[gps_coord$cc ==i & gps_coord$pnl_urban == 1], 
         col = 'red', pch = 20, cex = cex)
  points(gps_coord$lon[gps_coord$cc ==i & gps_coord$pnl_urban == 0], 
         gps_coord$lat[gps_coord$cc ==i & gps_coord$pnl_urban == 0], 
         col = 'darkgreen', pch = 20, cex = cex)
}
dev.off()



